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Abstract. Various methods of treating outer boundaries in numerical relativity 
arc compared using a simple test problem: a Schwarzschild black hole with 
an outgoing gravitational wave perturbation. Numerical solutions computed 
using different boundary treatments are compared to a 'reference' numerical 
solution obtained by placing the outer boundary at a very large radius. For 
each boundary treatment, the full solutions including constraint violations and 
extracted gravitational waves are compared to those of the reference solution, 
thereby assessing the reflections caused by the artificial boundary. These tests are 
based on a first-order generalized harmonic formulation of the Einstein equations 
and are implemented using a pseudo-spectral collocation method. Constraint- 
preserving boundary conditions for this system are reviewed, and an improved 
boundary condition on the gauge degrees of freedom is presented. Alternate 
boundary conditions evaluated here include freezing the incoming characteristic 
fields, Sommerfeld boundary conditions, and the constraint-preserving boundary 
conditions of Kreiss and Winicour. Rather different approaches to boundary 
treatments, such as sponge layers and spatial compactification, are also tested. 
Overall the best treatment found here combines boundary conditions that 
preserve the constraints, freeze the Newman- Penrose scalar 'to, and control gauge 
reflections. 



PACS numbers: 04.25.Dm, 02.60.Lj, 02.60.Cb 

1. Introduction 

A fundamental problem in numerical relativity is the need to solve Einstein's equations 
on spatially unbounded domains with finite computer resources. There are various 
ways of addressing this issue. Most often, the spatial domain is truncated at a finite 
distance and suitable boundary conditions are imposed at the artificial boundary. 
A different approach is to compactify the domain by using spatial coordinates that 
bring spatial infinity to a finite location on the computational grid. Another method 
often used for wave-like problems (although it is not commonly used in numerical 
relativity) includes so-called sponge layers which damp the waves near the outer 
boundary of the computational domain. The purpose of this paper is to compare these 
various methods by testing their ability to accurately reproduce dynamical solutions 
of Einstein's equations. 

An ideal boundary treatment would produce a solution to Einstein's equations 
that is identical (within the computational domain) to the corresponding solution 
obtained on an unbounded domain. In particular, no spurious gravitational radiation 
or constraint violations should enter the computational domain through the artificial 
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boundary. We can use this principle to test the various boundary treatments in the 
following way. First we compute a reference solution using a very large computational 
domain, large enough that its boundary remains out of causal contact with the interior 
spacetime region where comparisons are being made. Next we compute the same 
solution using a domain truncated at a smaller distance where one of the boundary 
treatments is used: we either impose boundary conditions there, compactify spatial 
infinity, or add a sponge layer. Finally we compare the solution on the smaller domain 
with the reference solution, measuring the reflections and constraint violations caused 
by the boundary treatment. Assessing boundary conditions by comparing with a 
reference solution on a much larger domain or a known analytic solution is a common 
practice in computational science. For applications to numerical relativity see e.g. [1], 
chapter 8 of [2] , and [3j HI [5] . 

The particular test problem used in this paper is a Schwarzschild black hole 
with an outgoing gravitational wave perturbation. The interior of the black hole is 
excised; all the characteristic fields propagate into the black hole (and out of the 
computational domain) at the inner boundary and hence no boundary conditions 
are needed there. Our numerical implementation uses a pseudo-spectral collocation 
method. See |Appendix A] for details on the initial data, the numerical methods, and 
the quantities that we compare between the solutions. 

We perform all of these tests using a first-order generalized harmonic formulation 
of the Einstein equations (see [5] and references therein) . In section [2] we discuss 
the construction of boundary conditions for this system that prevent the influx of 
constraint violations, and that limit the spurious incoming gravitational radiation 
by controlling the Newman-Penrose scalar ^ a t the boundary. We also improve 
the boundary conditions on the gauge degrees of freedom by studying small gauge 
perturbations of flat spacetime. We then evaluate the performance of these boundary 
conditions on our test problem: measuring the reflections and constraint violations 
caused by the computational boundary, and determining how these reflections vary 
with the radius of the boundary. 

Section [3] evaluates the performance of a variety of other widely used boundary 
conditions on our test problem. First we test the simple boundary conditions that 
freeze all the incoming characteristic fields at the boundary. We also test the 
commonly used variant of this, the Sommcrfcld boundary conditions, used in many 
binary black hole simulations [Tj [HI EH EH EE] based on the BSSN [TJ |TJ] formulation 
of Einstein's equations. Finally in section [3] we evaluate the constraint-preserving 
boundary conditions proposed by Kreiss and Winicour |14j . which differ from those 
discussed in section [2] mainly by our use of a physical boundary condition that controls 
•Jo- 
in section [4] we evaluate two boundary treatments that are alternatives to 
imposing local boundary conditions at a finite outer boundary. The first is the spatial 
compactification method used e.g. by Pretorius [T5J EH E2 m his groundbreaking 
binary black hole evolutions. In this treatment a coordinate transformation maps 
spatial infinity to a finite location on the computational grid. As waves travel out, 
they become increasingly blue-shifted with respect to the compactified coordinates 
and ultimately they fail to be resolved. Hence numerical dissipation is applied, which 
damps away these short-wavelength features. We measure the reflections and the 
constraint violations generated by the waves in our test problem as they interact with 
this boundary treatment. Finally in section [4] we implement and test a sponge layer 
method for Einstein's equations. 
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One of the main objectives of current binary black hole simulations is the 
computation of reliable waveforms for gravitational wave data analysis. Therefore it is 
important to evaluate how the various boundary treatments affect the accuracy of the 
extracted waveforms. In section|5] we compute the Newman-Penrose scalar ^4 (which 
describes the outgoing waves) on an extraction sphere close to the outer boundary (or 
compactified region, or sponge layer, respectively) and compare it with the analogous 
^4 from the reference solution. We also compare the measured reflections caused 
by our ty controlling boundary condition with the analytical predictions of these 
reflections made by Buchman and Sarbach [151 119j . 

Finally we discuss the implications of our results in section [6] and we also describe 
briefly a number of other boundary treatments which we do not test here. 

2. Constraint-preserving boundary conditions 

In this section, we briefly review the generalized harmonic form of the Einstein 
evolution system used in our tests. The method of constructing constraint-preserving 
boundary conditions (CPBCs) for this system is also discussed, and an improved 
boundary condition for the gauge degrees of freedom is derived. The numerical 
performance of these boundary conditions is evaluated using our test problem, and 
the dependence of the spurious reflections as a function of the boundary radius is 
measured. 

2.1. The generalized harmonic evolution system 

The formulation of Einstein's equations employed here uses generalized harmonic 
gauge conditions, in which the coordinates x a obey the wave equation 

Ox a = H a (x,ip), (1) 

where □ = ijj ab (d a d b — T c ab d c ) is the covariant scalar wave operator, with ip ab the 
spacetime metric and T c ab the associated metric connection. In this formulation of 
the Einstein system the gauge source function H a may be chosen freely as a function 
of the coordinates and of the spacetime metric ij) ab (but not derivatives of ip ab ). 

As is well known, the Einstein equations reduce to a set of coupled wave equations 
when the gauge is specified by equation (JlJ. We write this system in first-order form, 
both in time and space, by introducing the additional variables Q iab = di^Pab and 
n a b = —t c d c ip a b, where t c is the future directed unit normal to the t = const, 
hypersurfaces. Here lower-case Latin indices from the beginning of the alphabet 
denote four-dimensional spacetime quantities, whereas lower-case Latin indices from 
the middle of the alphabet are spatial. The principal parts of these evolution equations 
are given by [J] 

d t 4>ab ^ 0, 

d t Ii ab ~ N k d k IL ab - Ng kl d k ^ ab - l2 N k d k ^ ab , (2) 
d t $ iab ~ N k d k $ lab - Nd z Tl ab + N^O^ab, 

where ~ indicates that purely algebraic terms have been omitted, gij is the spatial 
metric of the t = const, slices, and N and N l are the lapse function and shift vector, 

X The parameter 71 of [6] is chosen to be —1, which ensures that the equations are linearly degenerate. 
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respectively. The parameter 72 was introduced in [6J in order to damp violations of 
the three-index constraint 

C iab = diijj ab - $ iab = 0. (3) 

We also include terms of lower derivative order that arc designed to damp violations 
of the harmonic gauge constraint [20] 

C a = -Dx a + H a = V b T a6c + H a = 0. (4) 

The system ([2| is symmetric hyperbolic. The characteristic fields in the direction 
Hi (where n a t a — 0) are given by 

u° ab = i/>ab, speed 0, (5) 

u£ = II o6 ± $„ a6 - 72^6, speed - N n ± N, (6) 

u 2 Aab = ®Aab, speed - N n . (7) 
For future reference, we also define 

u)± = n a6 ± $„ ab . (8) 

Here and in the following, an index n denotes contraction with Uj, while upper- 
case Latin indices A,B,... are orthogonal to n, e.g. va = Pa%v 1 where P a b = 
ip a b — n a rib + t a tb- For further details, we refer the interested reader to [B]. 

2.2. Construction of boundary conditions 

Our construction of boundary conditions for the generalized harmonic evolution 
system can be divided into three parts: constraint-preserving, physical, and gauge 
boundary conditions. 

In order to impose constraint-preserving boundary conditions, we derive the 
subsidiary evolution system that the constraints ^ and Q obey as a consequence of 
the main evolution equations ([2| . The incoming modes of the subsidiary system are 
then required to vanish at the boundary (cf. [211 E21 ESS Ell E3 ESI E3 EU ES]) . For 
instance, the harmonic gauge constraint Q obeys a wave equation 

□C a = (lower-order terms homogeneous in the constraints) (9) 

and the corresponding incoming fields will involve first derivatives of C a . In terms of 
the incoming modes it ao ~ (6| of the main evolution equations, the resulting constraint- 
preserving boundary conditions can be written in the form 

Pab^dnU^ = (\P ab P Cd - 2l (a P b) {c k d) + l a l b k C k d )d n ulj 

= (tangential derivatives), (10) 

where P c is a projection operator of rank 4 (cf. [BJ). Here nt now refers to the outward- 
pointing unit spatial normal to the boundary, l a = (t a + n a )/^/2, k a — (t a ~ n a )/^/2, 
and = denotes equality at the boundary. If the shift vector points towards the exterior 
at the boundary {N n >0), the fields u 2 Aab (JT]) are incoming as well and we obtain a 
boundary condition on them by requiring the components C n Aab of the four-index 
constraint 

C ijab = -2d^ j]ab (11) 

to vanish at the boundary. 

An acceptable physical boundary condition should require that no gravitational 
radiation enter the computational domain from the outside (except for backscatter 
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off the spacetime curvature, an effect that is a first-order correction in M/R). 
Gravitational radiation may be described by the evolution system that the Weyl tensor 
obeys by virtue of the Bianchi identities (see e.g. |27j). Our boundary condition 
requires the incoming characteristic fields of this system to vanish at the outer 
boundary. These incoming fields are proportional to the Newman-Penrose scalar vf^ 
(evaluated for a Newman-Penrose null tetrad containing the vectors l a and k a ). Hence 
the physical boundary condition we use is [271 E21 [30j [29l [31] 

9t*o = 0, (12) 



which can be written in a form similar to ( 10 1 



Pl b cd d n ul2 = (P a c P b d - \p ah p cd )d n u\r d 

= (tangential derivatives). (13) 
Here P p is a projection operator of rank 2 that is orthogonal to P c [BJ . We remark that 



( 12 I still causes some, albeit very small, spurious reflections of gravitational radiation. 
It can be viewed as the lowest level in a hierarchy of perfectly absorbing boundary 
conditions for linearized gravity [18| H"9]. 



The constraint-preserving (10) and physical (13) boundary conditions together 
constrain six components of the main incoming fields uzT. The remaining four 
components correspond to gauge degrees of freedom. In the past we chose simply 
to freeze those components in time [BJ, 

Pa b cd dtu^ = 0, (14) 

where P G = I-P C -P P . 

The initial-boundary value problem (IB VP) for the boundary conditions discussed 
so far was shown in |32j to be boundary-stable, which is a (rather strong) necessary 
condition for well posedness. These boundary conditions have been successfully used 
in long-term stable evolutions of single and binary black hole spacetimes [BJ [33 , 34J . In 
the following subsection, we present an improvement to the gauge boundary condition 



(14 1 motivated by the evolution of gauge perturbations about flat spacetime. 



2.3. Improved gauge boundary condition 

Let us assume that near the outer boundary, the spacetime is close to Minkowski space 
in standard coordinates (H a = 0) so that the Einstein equations may be linearized 
about that background. This assumption is reasonable because for the dominant 
wavenumber of the outgoing pulse (k = 1.6/M) and the boundary radius we typically 
consider (R = 41. 9M), we have kR 3> 1 and R ^> M. Furthermore, we assume that 
the outer boundary is a coordinate sphere of radius r = R. 

We begin by noting that harmonic gauge does not fix the coordinates completely: 
infinitesimal coordinate transformations 

x a ^x a + C (15) 
are still allowed provided the displacement vector satisfies the wave equation, 

DC = 0. (16) 
Under such a coordinate transformation, the metric changes by 

<tya6 = -2fl(„&). (17) 



A closer inspection [32] of the projection operator P G in (14 1 shows that the gauge 



boundary conditions control the components l a 5^jj a b of the perturbations, where 
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l a = (t a + ri a )/v2 is the outgoing null vector normal to the boundary. It is 
interesting to observe that these components vanish in the ingoing radiation gauge 
|35j . However, imposing radiation gauge on the entire spacetime is not possible in 
spacetimes containing strong-field regions, which will always generate perturbations 
l a Sipab that propagate into the far field. A reasonable condition to require then is that 
these perturbations pass through the boundary without causing strong reflections. 
Each Cartesian component of the vector l a 6ip a b obeys the scalar wave equation 

□V> = 0. (18) 

Solutions to this equation can be written in the form 



i 



^ = 2^ 2^ *WM)lM*,r), (19) 

l — l m— — l 

where the Y/ TO are the standard spherical harmonics and the ipi are linear combinations 
of outgoing (+) and incoming (— ) solutions 

^r)=r l - 1 (^j Ft{TTt), (20) 

Fjr(x) being arbitrary functions. A boundary condition is needed on ip that eliminates 
the incoming part of these solutions. In |36j . a hierarchy of boundary conditions 
is constructed that accomplish this task for all I ^ L. This idea was applied to 
the evolution of the Weyl curvature in |18j in order to construct improved physical 
boundary conditions. For the gauge boundary conditions considered here, we restrict 
ourselves to the L = member of the hierarchy, which corresponds to the Sommerfeld 
condition [|] 

(d t + d r + r- l )i) = Q. (21) 
In contrast, our old gauge boundary condition that froze the incoming characteristic 



field, as in (14), is given by 

(d t + d r + 72)^ = 0, (22) 
where 72 is the constraint damping parameter. 



This Sommerfeld boundary condition (21 1 is much less reflective than the freezing 



condition (22 1. To see this, we consider a solution of the form 

i>i = + pii>r (23) 

with generating functions 

Ff(x) = e ±ikx , (24) 

where k € K is the wave number. Substituting this solution into the boundary 
conditions (211 resp. (22 1, we solve for the reflection coefficient pi. Figure [I] shows 



\pi I for a typical range of wave numbers k and outer boundary radii R used for the 
numerical tests in this paper. (The dominant wave number of the outgoing pulse is 
k sw 1.6/M and in most cases, we place the outer boundary at R — 41.9M.) We 
see that \pi\ is much smaller (by about 3 orders of magnitude) for the Sommerfeld 
condition than for the freezing condition. 

§ To avoid confusion, we remark that in |5l ll4| , the term 'Sommerfeld condition' is used in reference 
to a condition of the form (dt + d r )u = 0, i.e. without the extra r~ 1 term due to our polar coordinates. 
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Figure 1. Predicted reflection coefficients p; for freezing (dotted) and 
Sommerfeld (solid) boundary conditions as functions of wave number k and outer 
boundary radius R. The curves for different I are visually indistinguishable in the 
freezing case. Note also that po = for the Sommerfeld condition. 



In the notation of the previous subsection, the improved gauge boundary 



condition (21 1 reads (after taking a time derivative), 



^ cd 9tfe" + (72-r-> cd ] = 0. 



(25) 



We remark that the extra terms in (25) as compared with the old condition (14 1 are 



of lower derivative order, so that the high-frequency stability result of 
immediately to these modified gauge boundary conditions. 



extends 



2.4- Numerical results 

The numerical tests of the various boundary conditions performed in this paper are 
described in some detail in Appendix A Figure [2] compares the numerical performance 
of our new CPBCs (flO), (fl~3j), p5] with our old ones (|l0|), (fTTj), (pj), (fllf . The 



outer boundary is placed at radius R — 41. 9M for these particular tests. Shown are 
the discrete L°° and L 2 norms of the difference AU between the numerical solution 



and the reference solution, and also the violations of the constraints C (see Appendix 
A. 4 for precise definitions of these quantities). The reference solution has an outer 
boundary at radius 961. 9M and is computed using our old CPBCs; thus for t < 920Af 
the outer boundary of the reference solution is out of causal contact with the region 
where AU and C are computed. 

In the difference AU we see a reflection that originates when the wave reaches the 
boundary at t ~ R and then amplifies as it moves inward in the spherical geometry, 
assuming its maximum at t ps 2R. This feature is much more prominent in the L°° 
norm than in the L 2 norm, which is why we display only the L°° norm in subsequent 
plots. The reflection is much smaller (by a factor of ~ R/M) for the new boundary 
conditions as compared with the old ones. Even at later times, the new boundary 
conditions result in a smaller AU, which in contrast to the old conditions appears to 
decrease as resolution is increased. 

We would like to stress that AU is a coordinate dependent quantity. Hence a 
smaller AU does not necessarily mean that the boundary treatment is 'better' in a 
physically meaningful sense. If however the aim is to produce a solution that is as 
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Figure 2. Old (solid) vs. new (dotted) CPBCs. Four different resolutions are 
shown: (N r ,L) = (21,8), (31,10), (41,12), and (51,14). The outer boundary is 
at _R = 41. 9M. 



close to the reference solution in the same coordinates, the choice of gauge boundary 
conditions does become important. Gauge reflections can in principle also impair the 
numerical accuracy of gauge-invariant quantities because much numerical resolution 
is wasted on resolving the gauge reflections. This is particularly the case when the 
gauge excitations in question are high-frequency modes such as those produced along 
with the so-called 'junk radiation' in binary black hole initial data. 

There is no discernible difference between the two sets of boundary conditions as 
far as constraint violations are concerned, which is what we expect because both of 
them are constraint-preserving. 

We close this section by investigating the dependence of the reflections on the 
radius of the outer boundary (figure [3]). The amplitude of the first peak in ||Ai/||oo 
decreases as the boundary is moved outward, roughly like l/R. At late times, there 
appears to be a power-law growth of that quantity at a rate that increases slightly 
with resolution. Inspection of the constraints (also in figure |3| and ^4 (figure 10 1 
suggests that this is a pure gauge effect. This blow-up is completely dominated by 
the innermost domain, which contains a long-wavelength feature that is growing in 
time. We speculate that this problem might be cured by a more clever choice of gauge 
source function close to the black hole horizon. 
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Figure 3. New CPBCs at different radii. Top half: all radii at the highest 
resolution, bottom half: R = 121. 9M at all resolutions. In the top right panel, 
curves for all outer boundary radii coincide. 



3. Alternate boundary conditions 

In this section, we consider several alternate boundary conditions that are often used 
in numerical relativity. All of these are local conditions imposed at a finite boundary 
radius, then in section [4] we consider some additional non-local boundary treatments. 
We run the alternate boundary conditions on our test problem and compare the results 



with the CPBCs (using the new gauge boundary condition (25 1 ) 



3.1. Freezing the incoming fields 

A very simple boundary condition is obtained by freezing in time all the incoming 
fields at the boundary, i.e., 

dtU 1 ^ = (and d t u 2 Aab = if N n > 0) . 



(26) 



This boundary condition is attractive from a mathematical point of view because it 
is of maximally dissipative type and hence, together with the symmetric hyperbolic 
evolution equations yields a strongly well-posed IB VP [37J [Ml ■ However, in 
general this boundary condition is not compatible with the constraints. 
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Figure 4. Freezing (solid) vs. new CPBCs (dotted). Four different resolutions 
are shown: (N r ,L) = (21,8), (31,10), (41,12), and (51,14). For freezing 
boundary conditions, both ||AW|| and C converge to a nonzero function with 
increasing resolution. The outer boundary is at R = 41.9Af. 



The left side of figure [4] demonstrates that freezing boundary conditions cause a 
significantly larger (by ss 3 orders of magnitude) initial reflection than our CPBCs. 
The difference with respect to the reference solution remains large in the subsequent 
evolution and unlike for the CPBCs does not decrease with increasing resolution. 
Furthermore, the violations of the constraints (right side of figure |4| do not converge 
away. This means that a solution to the Einstein equations is not obtained in the 
continuum limit. 



3.2. Sommerfeld boundary conditions 

A boundary condition that is often imposed in conjunction with the BSSN [T31 E3] 
formulation of the Einstein equations is a Sommerfeld condition on all the components 
of the spatial metric gij and extrinsic curvature Kij , 

(dt + dr + r- 1 ) ( 9ij ~^ =0. (27) 

This condition has been used for example in many recent binary black hole 



simulations [3 |8j |9j [10l [TT] . We cannot impose precisely the conditions ( 27 ) in our 
simulations because there is no one-to-one relationship between and Kij, and the 
incoming characteristic fields of our generalized harmonic formulation of Einstein's 
equations. Instead we consider the similar condition 

(d t + d r + r- 1 )^ - Vab ) = (28) 

on all the components of the spacetime metric {rj ab being the Minkowski metric). A 
very similar boundary condition (without the term) has recently been used in the 
generalized harmonic evolutions of [4"U] , 

In our formulation, boundary conditions are required not on the spacetime metric 
itself but only on certain combinations of its derivatives. By taking a time derivative 



of (28 1 and rewriting in terms of incoming characteristic fields, we obtain 

fitkb +(l2~r- 1 )i; ab }=0. (29) 




Figure 5. Sommerfeld (solid) vs. new CPBCs (dotted). Four different 
resolutions are shown: (N r ,L) = (21,8), (31,10), (41,12), and (51,14). The 
outer boundary is at R = 41. 9M. 



This then is our version of the Sommerfeld boundary condition (cf. (25 l), to be imposed 
on a spherical boundary in the far field (where linearized theory is assumed to be valid). 

Because the BSSN formulations using ( p7| are usually second-order in space, 
there is no analogue of our three-index constraint ([3| in that system. To mimic this 
situation in our tests of equation ( |29| ) , we also impose a CPBC on u\ ab as discussed in 
section |2.2| which together with our constraint damping terms ensures that violations 
of the three-index constraint (|3| are exponentially damped. 

Our version of Sommerfeld boundary conditions performs similarly on our test 
problem (figure [5]) to the freezing boundary conditions ( p6| (figure EJ. The initial pulse 
of reflections is smaller by w 2 orders of magnitude, but later ||AW|| grows to a similar 
level as for freezing boundary conditions. Again the constraints do not converge away, 
although this non-convergence appears only at somewhat higher resolutions than in 
the freezing case. 



3.3. Kreiss-Winicour boundary conditions 

Recently, Kreiss and Winicour [TJ] proposed a set of 'Sommerfeld- like' CPBCs for the 
harmonic Einstein equations and showed that they result in an IB VP that is well-posed 
in the generalized sense. Their boundary conditions were implemented and tested in 
[5]; here we compare their performance with the various other boundary treatments. 

The Kreiss-Winicour boundary conditions are obtained by requiring the harmonic 
constraint to vanish at the boundary, 

C a = 0. (30) 

In our notation, this can be written as an algebraic condition on part of the incoming 
fields u 1 ", 

P^' cd u^ = F a , (31) 

where 

p a c ' cd = ^[2fc( c <y*) -M> cd ], 

Fa = 4 l Kt - ^la^ bc ul+ + /"•'//;,., - \P a H hc ul c (32) 
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- l2t a + H a . 

The range of the projection operator P c is identical with that of P c defined in ( 10 ). 
For the unconstrained incoming fields u l ~ (i.e. u 1- without the 72 term, equation 
(JsJ ) , Kreiss and Winicour [2] specify certain free boundary data q^ b and q^ b . In our 
notation, 

pPcd-l- _ P pGcd~l- _ G 

r ab u cd ~ labi r ab U cd ~ lab- \ M ) 

In the linearized wave and gauge wave tests of [5], these boundary data are obtained 
from the known exact solutions. In the absence of an exact solution, it is suggested 
that the data could be obtained from an exterior Cauchy-characteristic or Cauchy- 
perturbative code. However, since we do not have such an exterior code, we compute 
the boundary data from the background solution, i.e. Schwarzschild spacetime. As in 
the Sommcrfeld case (section 3.2), we use a constraint-preserving boundary condition 
on u\ ab to emulate the second-order formulation of (HE], and this value of u\ ab is 
then used to compute F a in (32 1. 

Figure [6] shows the numerical results for our test problem. The magnitude of the 
initial reflections lies between that of freezing and Sommerfeld boundary conditions 
and is somewhat smaller at later times, though still larger than for our CPBCs at the 
higher resolutions. The constraints converge away with increasing resolution, as they 
should for a boundary condition that is consistent with the constraints. In a numerical 
simulation, violations of the constraints are in general present in the interior of the 
computational domain. These propagate as described by the constraint evolution 
system ^ and some may hit the outer boundary. The Dirichlet boundary conditions 
d30| might be expected to cause more reflections of constraint violations than our 



no-incoming-field conditions (10 1, however, no indications of this are seen in figure 
[6] Probably the constraint damping we use is sufficiently effective in eliminating the 
source of these reflections. 

We shall see in section 5.1 that the Kreiss- Winicour boundary conditions also 
cause larger errors in the physical degrees of freedom than our CPBCs. Since the 
main difference between the two sets of boundary conditions is our use of a physical 
boundary condition d t ^o = 0; we conclude that such a condition is crucial in reducing 
the reflections from the outer boundary. 



4. Alternate approaches 

So far we have only considered boundary conditions that are local algebraic or 
differential conditions imposed at the boundary of some finite computational domain. 
There are of course many ways of treating the outer boundary that do not fall into that 
category. In this section, we evaluate two such approaches: spatial compactification 
and sponge layers. 



4-1- Spatial compactification 

Spatial compactification is a method that has been widely used in numerical relativity, 
for instance in |41l 142] or more recently in the generalized harmonic binary black hole 
simulations of Pretorius [THl EE HZ] . 

The basic idea is to introduce spatial coordinates that map spacelike infinity to a 
finite location. Here we consider mappings that are functions of coordinate radius only 
(whereas Pretorius applies the mapping to each Cartesian coordinate separately) . We 




Figure 6. Kreiss-Winicour (solid) vs. new CPBCs (dotted). Four different 
resolutions are shown: (N r ,L) = (21,8), (31,10), (41,12), and (51,14). The 
outer boundary is at R = 41. 9M. 



have used two such mappings, named Tan and Inverse, as detailed in |Appcndix B.l| 
Each map has a scale R across which the mapping is (essentially) linear. The outermost 
grid point is placed at a very large but finite uncompactified radius (r = 10 17 M). With 
respect to the compactified radial coordinate, the characteristic speeds are below 
numerical roundoff there and hence no boundary condition should be needed. The 
following results were produced using constraint-preserving boundary conditions; we 
have checked for one simulation that using no boundary condition at all yields results 
that are visually indistinguishable from the ones presented here on the scales of figures 
[7][8] and[lO| 

As the waves travel outward, they become more and more blue-shifted with 
respect to the computational grid and are eventually no longer properly resolved. 
However, some form of artificial numerical dissipation is applied that acts as a low- 
pass filter and causes the waves to be damped as they become increasingly distorted. 
We have experimented with various such filters; see |Appcndix B.l for details. One of 
them (referred to as number 2 in the following) is designed to emulate as closely as 
possible the fourth-order Kreiss-Oliger dissipation used by Pretorius. 

In the following numerical comparisons, we evaluate the differences with respect 
to the reference solution only in the part of the domain where the compatification map 
is essentially linear, i.e. for r ^ R. First we compare the various filtering methods at 
a fixed resolution, using the Tan compactification mapping (figure [7j. The filters that 
are applied to the right side of the evolution equations (numbers 1 and 3, cf. table Bl I 
do somewhat better than those applied to the solution itself (numbers 2 and 4), and 
the Exponential filters (numbers 3 and 4) are slightly better than the Kreiss-Oliger 
filters (numbers 1 and 2). All of them are outperformed by the CPBCs (imposed at 
r = R). For our closest approximation to the dissipation used by Pretorius (number 
2), 1 1 AW 1 1 is comparable to constraint-preserving boundary conditions at the peak of 
reflections (at t 2R) but becomes larger by about 2 orders of magnitude at later 
times. The compactification methods also generate considerable constraint violations. 

Next we focus on the best filter (number 4) of the previous test but vary the 
resolution (figure |8j. We do see convergence of ||At/|| initially but the convergence 
degrades at later times. This is surprising at first because with increasing resolution, 
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Figure 7. TAN compactification with various filters vs. new CPBCs. Only the 
highest resolution (N r ,L) = (51,14) is shown. The compactification scale (and 
the radius of the outer boundary in the CPBC case) is R = 41. 9M. 
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Figure 8. TAN compactification with filter 4 (solid) vs. new CPBCs (dotted). 
Four different resolutions are shown: (N r ,L) = (21,8), (31,10), (41,12), and 
(51, 14). The compactification scale (and the radius of the outer boundary in the 
CPBC case) is R = A1.9M. 



the waves travel a longer distance before they fail to be resolved. Note however that 
the high-frequency filter is applied at each time step, as is done in the simulations 
of Pretorius. For higher resolutions, the time steps are smaller because of the CFL 
condition and the filter is applied more often, thus leading to a stronger damping of 
the waves. This may well lead to the observed loss of convergence with increasing 
resolution. The constraints appear to converge away in this test, although from figure 
[8] it appears that this will not persist for even higher resolutions. 

We have also evaluated the Inverse mapping described in |Appcndix B.l| The 
results are similar, but somewhat worse than the Tan mapping results shown here. 
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4.2. Sponge layers 

A method that has been used for a long time in computational science, in particular 
for spectral methods (see e.g. section 17.2.3 of [12] and references therein), involves 
so-called sponge layers. A sponge layer is introduced by modifying the evolution 
equations according to 

d t u= ...-7(r)(u-uo), (34) 

where uq refers to the unperturbed background solution (Schwarzschild spacetime in 
our case) and the smooth sponge function 7(7") > is large only close to the outer 
boundary of the computational domain. (Here we use uncompactified coordinates as 
in sections[2]and[3]) In this way, the waves are damped exponentially as they approach 
the outer boundary. Details on our particular choice of 7(7-) can be found in | Appendix 

E2 

We compare the sponge layer method with our CPBCs in figure [9j For the 
CPBCs, the boundary is either placed at R = 41. 9 M (the outer edge of the sponge- 
free region) or at R = 121.9M (the outer edge of the sponge). At early times (t < 2R), 
the 1 1 AW I |oo of the sponge layer method lies between that of the CPBCs for the two 
choices of outer boundary radius, whereas at later times, it is much larger than both 
versions of CPBCs. The constraint violations in the sponge runs do not converge 
away. 



5. Physical gravitational waves 

Perhaps the most important predictions of numerical relativity simulations at the 
present time are the gravitational waveforms produced by astrophysical systems like 
binary black holes. It is important therefore to understand how the accuracy of these 
waveforms is affected by the choice of boundary treatment. Physical gravitational 
radiation can be described by the Newman-Penrose scalars 'J 4 and ^o- The scalar 'J 4 
is dominated by the outgoing radiation (its ingoing part is suppressed by a factor of 
(fcr) 4 , where k is the wavenumber), whereas "4>q is dominated by the ingoing radiation 
(its outgoing part is suppressed by a factor of (fcr) 4 ). In this section we compare 
the gravitational waves extracted from the various boundary treatment solutions on 
a sphere of radius r = i? cx , using the methods described in | Appendix A. 5) 

We note that ^4 (^0) has a coordinate-invariant meaning only in the limit as 
future (past) null infinity is approached. The quantities computed at finite radius r will 
differ from those observed at infinity by terms of the order 0(l/r). In the particular 
case of perturbed Schwarzschild spacetime considered here, a gauge-invariant wave 
extraction method does exist even at finite radius (see e.g. [H] and references therein) 
but we do not adopt it here. Our purpose in this paper is merely to measure the 
effects on ^4 caused by the various boundary treatments. 



5.1. Difference of ^4 with respect to the reference solution 

We begin by evaluating A^ 4 = ^4 — v I / | cf , where 'J 4 is the Newman-Penrose scalar 
computed using one of the various boundary methods and v?5~ is the same quantity 
computed from the reference solution at the same extraction radius. The curves shown 
in figure 10 plot the maximum value of IAW4I over time intervals of length 20M (this 
time filtering averages over the high frequency quasi-normal oscillations of the black 
hole), normalized by the maximum value of |A^| over the entire evolution. The 
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Figure 9. Sponge layer method (solid) vs. new CPBCs at two different radii 
(dotted). Four different resolutions are shown: (jV r , L) = (21, 8), (31, 10), (41, 12), 
and (51, 14). The size of the sponge-free region is R = 41. 9M and |[AW||oo is only 
computed for r sC R. 



radius of the outer boundary (or the compactification scale, or the size of the sponge- 
free region, respectively) used for these comparisons is R — 41. 9M, and the radiation 
is extracted nearby at R ex = 40M. 

The first peak in |A^4| seen in figure 10 arises as the wave in our test problem 
passes outward through the extraction sphere at t as R cx . This peak is caused by 
a presently unknown (probably gauge) interaction between the outer boundary (or 
compactified region etc.) and the spacetime near the extraction sphere. We have 
verified that this interaction and its influence on the peak in A ^4 goes away if we 
move the outer boundary (or the extraction surface) so that they are not in causal 
contact as the outgoing wave pulse passes the extraction surface. 

Some of the outgoing radiation is reflected off the boundary. Most of this reflected 
radiation is subsequently absorbed by the black hole, but some of it excites the hole, 
which then emits quasi-normal mode radiation of exponentially decaying amplitude. 
This exponential decay can be clearly seen for most of the boundary treatments. 

In the case of freezing boundary conditions, nearly all of the outgoing quasi- 
normal mode radiation is reflected from the boundary because the reflection coefficient 
is nearly 1 for the wave number of the dominant mode, k — 0.376/A/ (cf. figure [l]). 
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It then re-excites the black hole, which again radiates and so forth. On average the 
amplitude of the reflections remains roughly constant in time for this case. This 
behaviour is consistent with the result shown in figure 3 of [6] for a similar perturbed 
black hole simulation. 

For the Sommerfeld and Kreiss-Winicour boundary conditions, the reflections are 
much smaller but still considerably larger (by 2 to 3 orders of magnitude) than for 
our CPBCs. We attribute this difference largely to our use of the physical boundary 



condition ( 12 1 



The spatial compactification method has the largest difference |A^4|, particularly 
at early times t ~ R (about 4 orders of magnitude larger than for the CPBCs). We 
suspect that this may be a consequence of the use of artificial dissipation, as discussed 
in section |4~T1 

The sponge layer method has the smallest errors at early times. This is not 
surprising because the outer boundary of the sponge layer is much further out at 
R = 121. 9M. However at later times when the waves begin to interact with the 
sponge layer, this method causes reflections comparable in amplitude to those using 
Sommerfeld boundary conditions. 

We also note that at late times the level of (A^l decreases significantly with 
resolution for the CPBCs, but not generally for the other boundary treatments. 

We think it is remarkable that the maximum relative error in the extracted 
physical radiation is quite small (10~ 5 to 10~ 3 ) in these tests, even for the less 
sophisticated boundary treatments such as the freezing or Sommerfeld boundary 
conditions. This success is due in part to the fact that the extraction radius, 
R cx = 40M, for this test problem is about ten wavelengths (of the initial radiation 
pulse) away from the central black hole. Our results are likely to be more accurate 
than those from typical binary black hole simulations, which place the outer boundary 
at two or three wavelengths. This suggests that current binary black hole codes 
using, for instance, Sommerfeld boundary conditions, can still produce waveforms 
that are useful for some aspects of gravitational wave data analysis provided the 
outer boundary is placed sufficiently far out. Data analysis applications needing 
high precision waveforms, however, such as source parameter measurement or high- 
amplitude supermassive binary black hole signal subtraction for LISA, will need to 
use a more sophisticated boundary treatment that produces smaller errors in ^4. 

5.2. Comparison with the predicted reflection coefficient 

Buchman and Sarbach [T51 [T5] have recently developed a hierarchy of increasingly 
absorbing physical boundary conditions for the Einstein equations by analyzing 
the equations describing the evolution of the Weyl curvature on both a flat and 
a Schwarzschild background spacetime. Their analysis predicts, in particular, the 
reflection coefficient p (defined as the ratio of the ingoing to the outgoing parts of the 
solution) that arises from the d{^o = physical boundary condition that we use. 

For quadrupolar radiation (as in our numerical tests), this reflection coefficient is 
given by equation (89) of [TH], 

P (kR) = l{kR)- A + 0(kR)~ 5 , (35) 

where k is the wave number of the gravitational radiation and R is the boundary 



radius. (As explained at the beginning of section 2.3 we assume the background 



spacetime to be flat; effects due to the backscattering would only enter at 0(M/R).) 
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Figure 10. Difference of \E'4 for the various alternate methods (solid) vs. the 
new CPBCs (dotted). Two resolutions are shown: (N r ,L) = (31, 10) and (51, 14). 
The radius of the outer boundary (or the compactification scale, or the size of 
the sponge-free region, respectively) is R = A1.9M and the waves are extracted 
at _R ex = 40 M. 
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By evaluating Wo and W4 at the extraction radius of our test, we find that the ratio 
^0/^4 agrees with their predicted p to leading order in l/(kR). We note that the 
tetrad we use for wave extraction (Appendix A. 5) does not agree exactly with that 
of |18j . However, the tetrads do agree for the unperturbed Schwarzschild solution, 
so that the errors introduced into * an d W4 due to our different choice of tetrad 
are second-order small in perturbation theory and hence the comparison with |18j is 
consistent. 

For a numerical solution using our new CPBCs, we evaluate the Newman- Penrose 
scalars ^^(t) and &4,(t) on extraction spheres located 1.9M inside the outer boundary. 
In figure [TT] we plot the time Fourier transforms of these quantities. We also plot 
|(fci?)~ 4| I'4, which by the above argument should agree with W to leading order in 



l/(kR). Figure 11 shows that the numerical agreement is reasonably good: roughly 
at the expected level of accuracy. The overall dependence of the predicted reflection 
coefficient p on k and R is captured very well. We surmise that the levelling off of 
our numerical Wo for fc > 3 is due to numerical roundoff effects. (Note the magnitude 
of W a t those frequencies.) For radii R > 200M, W is at the roundoff level for all 
frequencies. 



6. Discussion 



The purpose of this paper is to compare various methods of treating the outer 
boundary of the computational domain. We evaluate the performance of several often- 
used boundary treatments in numerical relativity by measuring the amount of spurious 
reflections and constraint violations they generate. To this end, we consider as a test 
problem an outgoing gravitational wave superimposed on a Schwarzschild black hole 
spacetime. First we compute this numerical solution on a reference domain, large 
enough that the influence of the outer boundary can be neglected. Then we repeat the 
evolution on smaller domains using one of the boundary treatments, either imposing 
local boundary conditions, compactifying the domain using a radial coordinate map, 
or installing a sponge layer. We use a first-order generalized harmonic formulation 
of the Einstein equations, although these boundary methods can be applied to other 
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formulations as well. Wc believe our results arc fairly independent of the particular 
formulation used. 

Our main conclusion is that our version of constraint-preserving boundary 
conditions performs better than any of the alternate treatments that we tested. Our 
boundary conditions include a limitation on the influx of spurious gravitational waves 
by freezing the Newman-Penrose scalar vffp at the boundary. We also introduce and 
test an improved boundary condition for the gauge degrees of freedom. 

For some of the simple boundary conditions, such as freezing or Sommcrfcld 
conditions, we find constraint violations that do not converge away with increasing 
resolution. The continuum limit does not satisfy Einstein's equations in these cases. 
Most of the alternate boundary conditions also generate considerable reflections as 
measured by AU, the norm of the difference with respect to the reference solution. In 
many cases, these reflections do not decrease significantly with increasing resolution. 

The difference norm AU that we use to measure boundary reflections includes the 
entire spacetime metric, not just the physical degrees of freedom. It is important then 
to evaluate separately the effects of the various boundary treatments on the physical 
degrees of freedom. We use the extracted outgoing radiation as approximated by 
the Newman-Penrose scalar 'J 4 for this purpose. Here our conclusions are somewhat 
different. Rather surprisingly, most of the boundary methods we consider generate 
relatively small errors in ^4. This suggests that if gravitational waveforms are only 
needed to an accuracy of, say, 1% (which is comparable to the discrepancies between 
recent binary black hole simulations [45 ) then even the simple Sommcrfcld conditions 
might be good enough. (For those, we find relative errors ~ 10~ 5 .) The largest 
relative errors in ^4 we find (~ 10~ 2 ) occur with our implementation of the spatial 
compactification method used by Pretorius (X5J ESI EZ] ■ We attribute these largely to 
the use of artificial dissipation. Undesirable effects of dissipation might be somewhat 
less severe in binary black hole evolutions, which have much larger wavelengths 
(A ~ 20 — 100M) than ours (A ~ 4M). Our tests suggest that the errors in ^4 can 
be made to decrease significantly with resolution only by using more sophisticated 
constraint preserving and physical boundary conditions. The importance of using a 
physical boundary condition on "to is illustrated in particular by the difference between 
the performance of our boundary conditions and those of Kreiss and Winicour |14j . 

Some caveats regarding the interpretation of our results must be stated. First, 
the ratio of the dominant wavelength to the radius of the outer boundary is typically 
much larger for binary black hole evolutions (where X/R > 0.5) than for the simple 
test problem considered here (where X/R ~ 0.1). Boundary treatments generally work 
better for smaller X/R, i.e. when the boundary is well out in the wave zone. Hence the 
results presented here are likely to be more accurate than those from typical binary 
black hole simulations. Second, we use spectral methods rather than finite-difference 
methods, which are more commonly used in numerical relativity at this time. This 
complicates the implementation of the kind of numerical dissipation that is crucial for 
the spatial compactification method to work. While we have attempted to construct 
a filter that mimics the finite-difference dissipation as closely as possible, a direct 
comparison is clearly impossible. In finite-difference methods, the error introduced by 
the type of numerical dissipation considered here is below the truncation error. Hence 
tests similar to ours but performed with a finite-difference method would not be able 
to detect the effect of dissipation. 

There are several directions in which the present work could be extended. For 
large values of the outer boundary radius, we observe a non-convergent power- 
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law growth of the error in our test problem when constraint-preserving boundary 
conditions are used; the origin of this growth should be investigated further. It would 
be interesting to implement and test the hierarchy of physical boundary conditions that 
are perfectly absorbing for linearized gravity (including leading-order corrections due 
to the curvature and backscatter) found recently by Buchman and Sarbach [TBI US] ■ 
Our boundary conditions could also be tested using known exact solutions such as 
gauge waves, and comparisons could be made with the results found in [5j. 

For completeness we also mention a number of additional outer boundary 
approaches that were not addressed in this paper, but would also be interesting future 
extensions of this research. In [JHl U7J, boundary conditions for the full nonlinear 
Einstein equations on a finite domain are obtained by matching to exact solutions of 
the linearized field equations at the boundary. Alternatively, the interior code could 
be matched to an 'outer module' that solves the linearized field equations numerically 
[15J IIHI [SHI [5T] . Other approaches involve matching the interior nonlinear Cauchy code 
to an outer characteristic code (see |52j for a review) or using hyperboloidal spacetime 
slices that can be compactified towards null infinity (see [53] for a review). 

Appendix A. Details on the numerical test problem 

Appendix A.l. Initial data 

The initial data used for our numerical tests are the same as in |27j . The background 
solution is a Schwarzschild black hole in Kerr-Schild coordinates, 

2M 

ds 2 = -dt 2 + (dt + dr) 2 + dr 2 + r 2 dn 2 . (A.l) 

r 

Throughout the paper, M refers to the bare black hole mass of the unperturbed 
background. We superpose an odd-parity outgoing quadrupolar wave perturbation 
constructed using Tcukolsky's method [S3]. Its generating function is taken to be a 
Gaussian G(r) = Aexp[-(r - r ) 2 /w 2 ] with i = 4x 1CT 3 , r = 5M, and w = 1.5M. 
The full non-linear initial value equations in the conformal thin sandwich formulation 
are then solved to obtain initial data that satisfy the constraints [55] ■ This yields 
initial values for the spatial metric, extrinsic curvature, lapse function, and shift 
vector. We note that after the superposition, the resulting solution is still nearly 
but not completely outgoing. 

Our generalized harmonic formulation of Einstein's equations requires initial data 
for the full spacetime metric and its first time derivative. These can be computed from 
the 3+1 quantities obtained above, provided we also choose initial values for the time 
derivatives of the lapse function and shift vector. These initial time derivatives are 
freely specifiable and are equivalent to the initial choice of the gauge source function 
H a ; we choose d t N = and d t N l = at t = 0. 

Appendix A. 2. Numerical method 

We use a pseudospectral collocation method as described for example in [27] . 

The computational domain for the test problem considered here is taken to 
be a spherical shell extending from r = 1.9M (just inside the horizon) out to 
some r = R. This domain is subdivided into spherical-shell subdomains of extent 
A7' = lOAf. On each subdomain, the numerical solution is expanded in Chebyshev 
polynomials in the radial direction and in spherical harmonics in the angular directions 
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(where each Cartesian tensor component is expanded in the standard scalar spherical 
harmonics). Typical resolutions are N r G {21,31,41,51} coefficients per subdomain 
for the Chebyshev series and I ^ L with L £ {8, 10, 12, 14} for the spherical harmonics. 

We change the outer boundary radius R by changing the number of subdomains 
while keeping the width Ar of each subdomain fixed; this facilitates direct comparisons 
between runs with different values of R. For example, the innermost four subdomains 
of the reference solution (which has a total of 96 subdomains and R = 961. 9M) are 
identical to the four subdomains used to compute the solution with R = 41. 9M. 

The evolution equations are integrated in time using a fourth-order Runge-Kutta 
scheme, with a Courant factor At/Ax m ; n of at most 2.25, where Ax m j n is the smallest 
distance between two neighbouring collocation points. As described in [27], the top 
four coefficients in the tensor spherical harmonic expansion of each of our evolved 
quantities is set to zero after each time step; this eliminates an instability associated 
with the inconsistent mixing of tensor spherical harmonics in our approach. 

We use two methods of numerically implementing boundary conditions; the choice 
of method depends on the type of boundary conditions. Boundary conditions that can 
be expressed as algebraic relations involving the characteristic fields are implemented 
using a penalty method (see [SS] and references therein; in the context of finite- 
difference methods see also |57| and references therein). In particular, we use a penalty 



method to implement the Kreiss-Winicour boundary conditions (cf. section 3.3 1 
and to impose boundary conditions at the internal boundaries between neighbouring 
subdomains. Boundary conditions that are expressed in terms of the time derivatives 
of the characteristic fields are implemented using the method of Bj0rhus [58] , where 
the time derivatives of the incoming characteristic fields are replaced at the boundary 
with the relevant boundary condition. All boundary conditions in this paper besides 
those mentioned above are implemented using the Bj0rhus method. 



Appendix A. 3. Gauge source functions 

Our generalized harmonic formulation [6 j of Einstein's equations allows for gauge 
source functions that depend arbitrarily on the coordinates and the spacetime metric: 
H a = H a {t,x,ip). The generalized harmonic evolution equations are equivalent to 
Einstein's equations only if the constraint Q remains satisfied. 

We choose the time derivatives of lapse and shift to be zero at the beginning of 
the simulation; this determines the initial value of H a via the constraint Q . For the 
subsequent evolution, we hold this H a fixed in time. 



Appendix A. 4-- Error quantities 

We use two different measurements of the errors in our solutions, which we monitor 
during our numerical evolutions. First, given a numerical solution (ip a b, n a b, &i a b)i 
the difference between that solution and the reference solution (ip^ b , , Q^P) ^ s 
computed with the following norm at each point in space, 

AU = [6 ab S cd (M~ 2 Aip ac Aip bd + An QC An M 

+g^A^ ac A^ jbd )] 1/ \ (A.2) 

where Aip a i, means V>ab — i^ab^i an d similarly for AII a f, and A$ ia b. Second, we define 
a quantity C that measures the violations in all of the constraints of our system, 

C = [S^^a^b + g^iCiaCjb + g kl S cd CikacCjibd) 
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+M- 2 (C a C b + g ij 6° d C iac C jbd ))] 1/2 , (A.3) 

where J- a and Ci a are first derivatives of C a defined in [6] . To compute global error 
measures, a spatial norm 1 1 • 1 1, either the L°° norm or the L 2 norm, is applied separately 
to AU and C. 

The question often arises as to the significance of particular values of ||AW|| and 
||C||. For example, is a simulation with ||C|| = 10~ 2 good to one percent accuracy? To 
make it easier to answer such questions, we normalize both ||AW|| and ||C|| as follows, 
and we always plot normalized quantities. 

We divide ||AW|| by a normalization factor ||AWo||, defined as the difference 
between a given solution at t = and the unperturbed Schwarzschild background; 



i.e., the quantity ||AWo|| is computed from (A. 2) using the unperturbed Schwarzschild 



solution instead of the reference solution. Since ||AWo|| is evaluated at t = 0, it depends 
only on the initial data used in the simulation, and is a measure of the amplitude 
of the superposed gravitational wave perturbation. For the initial data used here, 
HAWolloo = 6 x 1CT 3 and ||AW || 2 = 1.4 x 1CT 4 . The quantity ||AW||/||AW || is more 
easily interpreted than ||At/||; for example, ||AW||/||AWo|| is unity when the difference 
from the reference solution is of the same size as the initial perturbation. 

Similarly, the constraint energy norm ||C|| is divided by the norm of the first 
derivatives \\dU\\ (at the respective time), 

du = [^s ab s cd (Ar 2 d^ ac d^ bd + ^n oc ^n M 

+g kl d i $ kac d j $i bd )] 1/2 . (A.4) 

The constraints for our system are linear combinations of the first derivatives of the 
fields, hence ||C||/||3W|| ~ 1 corresponds to a complete violation of the constraints. 

Appendix A. 5. Wave extraction 

For evaluating gravitational waveforms, we compute the Newman-Penrose scalars 

*o = -C abcd l a m b l c m d , * 4 = -C abcd k a m b k c m d , (A.5) 

where C abc d is the Weyl tensor, l a and k a are outgoing and ingoing null vectors 
normalized according to l a k a = — 1, m a is a complex unit null spatial vector orthogonal 
to l a and k a , and fh a is the complex conjugate of m a . For perturbations of flat 
spacetime, there is a standard choice for the vectors l a , k a , and m a . In general curved 
spacetimes, however, no such prescription for the tetrad exists that would produce 
coordinate-independent quantities and ^4 at finite radius. We choose the null 
vectors according to 

l a = i (t a + n a ), k a = (t a - n a ) , (A.6) 

where t a is the future-pointing unit timelike normal to the t = const, slices and n a is 
the unit spacelike normal to the extraction sphere. Finally, we choose 

1/9 .15 



V2r \d0 ' "sinC u 1 ' (AJ) 

where (r, 8, <f>) are spherical coordinates on the r — R cx — const, extraction sphere. 
Note that our choice of m a is not exactly null nor of unit magnitude at finite extraction 
radius. However, the tetrad is orthonormal for the unperturbed Schwarzschild 
solution, so that the errors introduced into and ^4 because of the lack of tetrad 
orthonormality will be second-order small in perturbation theory. 
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The quantity ^4 corresponds to outgoing radiation in the limit of r — > 00, 
t — r — const., i.e. as future null infinity is approached. Similarly corresponds 
to ingoing radiation as past null infinity is approached. At finite extraction radius, 
^4 and "fo will disagree with the waveforms observed at infinity by terms of the order 

O(iiex)- 1 . 

We decompose the quantities ^4 and Vto m terms of spin-weighted spherical 
harmonics of spin-weight —2 on the extraction surface. Since our perturbation is 
an odd-parity quadrupole wave, the imaginary part of the (I = 2, m = 0) spherical 
harmonic is by far the dominant contribution to ^4, and we only display that mode 
in our plots. We normalize the curves in our graphs by the maximum (in time) value 
of I $4 1 at the extraction radius i? eX) which for 7? ox = 40M is max = 6 x 1CP 4 . 



Appendix B. Details of the alternate approaches 

In this appendix, we provide some more details on the alternate boundary treatments 
discussed in section [3] spatial compactification and sponge layers. 



Appendix B.l. Spatial compactification 

We implement spatial compactification by introducing a radial coordinate 
transformation x — > r(x) that maps a compact ball on the computational grid with 
x € [0,x max ] to the full unbounded physical slice with r € [0,oo]. We consider two 
such mappings. The Tan mapping is similar to the one used by Pretorius (131 [Trjl [T7] 
and is given by 

—J, 0^x<2R. (B.l) 
The scale R determines the range in physical radius r across which the map is 



essentially linear (see figure Bll. When comparing compactification with other 
boundary treatments, we compare quantities only in the region r < R. (The scale 
R is equal to unity in the work of Pretorius. He uses mesh refinement to obtain the 
appropriate resolution close to the origin, while we fix the resolution and choose the 
scale R appropriately.) We also tested an Inverse map defined by 

r x, < x sC R , 

Inverse^) = { R 2 „ rtn ( B - 2 ) 
, R< x <2R, 

see figure Bl This map is only C 1 at x = R, but we maintain spectral accuracy in 
our tests by placing this surface at the boundary between spectral subdomains. 

Dissipation is needed to remove the short wavelength components of the 
waves as they travel outward on the compactified computational grid and become 
unresolved. We apply this dissipation only in the radial direction, but everywhere 
in the computational domain. In spectral methods, dissipation can be conveniently 
implemented in the form of a spectral filter. This filter is applied by multiplying each 



spectral expansion coefficient of index & by a function f(k). (See Appendix A. 2 for 
details on the pseudospectral method we use.) Higher values of k correspond to shorter 
wavelengths in the numerical approximation; let fc max be the highest index used in the 
spectral expansion. The first filter function we consider is the closest analogue in the 
context of our spectral methods to Kreiss-Oliger [53] dissipation, 

/kreiss-OligerW = 1 - e sin4 ( ^ — ) . < e < 1. (B.3) 
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or 
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Figure Bl. Compactification mappings (left) and filter functions (right). The 
dashed line indicates the boundary of the region in where the compactification 
mapping is (essentially) linear. 



Typical values of the parameter e used by Pretorius are e € [0.2, 0.5]; we use e = 0.25. 

This filter was derived via a comparison with finite-difference methods as follows. 
In the finite-difference approach, a numerical solution u is represented on a set of 
equidistant grid points Xj. (It suffices to consider the one-dimensional case here.) 
Some form of numerical dissipation is usually required for the finite-difference method 
to be stable. The one that is most often used for second-order accurate methods is 
fourth-order Kreiss-Oliger dissipation [59]. One possible implementation of this, used 
e.g. by Pretorius, amounts to replacing 

^ 4 £ 4 ) u (B.4) 

at each time step, where h is the grid spacing and D A is the second-order accurate 
centred finite difference operator approximating the fourth derivative, 

D A Ui = h~ A (uj-2 — 4uj-i + 6uj — 4uj + i + Wj+2)- (B-5) 

Taking u to be a Fourier mode Uj = exp(ikxj), it follows that the mode is damped 
by a frequency-dependent factor, 



Flu] 



1 



,(*0 



F[u 



(fc)l 



4 / 7rfc 



(B.6) 



where fc r 



n/(2h) is the Nyquist frequency. Thus we obtain the filter function 



(B.3I. 



Strictly speaking, the above analysis only applies to Fourier expansions and 
not to the Chebyshev expansions we use. Nevertheless, we apply the filter in the 
form (B.3l to our Chebyshev expansion coefficients. Note that in (B.6), each spectral 
coefficient is filtered separately; this is not true for the analogous calculation for 
a Chebyshev expansion. 

We also use a different filter function, which we call the Exponential filter, that 
is often used in spectral methods (see [BD] and references therein), 

k 



/e> 



XPONENTIAL 



(fc) = exp 



(B.7) 
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No. 


Type 


Parameters 


Applied to 


1 


Kreiss-Oliger 


e = 0.25 




right side 


2 


Kreiss-Oliger 


e = 0.25 




solution 


3 


Exponential 


(7 = 0.76, p = 


13 


right side 


4 


Exponential 


(7 = 0.76, p = 


13 


solution 



Table Bl. Details of the filtering methods 



Typical values of the parameters are a — 0.76 and p = 13. This choice of parameters 
gives less dissipation at small values of k than the Kreiss-Oliger filter, and also ensures 
that /(fcmax) ~ 10~ 16 is at the level of the numerical roundoff error. 

There are various ways the filters can be applied in a numerical evolution. We 
have experimented with two different methods. In the first method, the filter is applied 
to the right side of the equations, i.e. the evolution equations dtu — S are modified 
according to dtu — F[S], where F[S] is the filtered right side. In the second method, 
the filter is instead applied to the solution itself, i.e. after each substep of the time 



integrator (cf. Appendix A.2|, the numerical solution u is replaced with its filtered 
version F[u]. This second method is closest to how the Kreiss-Oliger filter is applied 
by Pretorius. 

For our numerical tests, we have used four different combinations of the various 
options described above. They are summarized in table |B1| 

Appendix B.2. Sponge layers 



For sponge layers we must specify a sponge profile function 7(7") , as defined in ( 34 1 . We 
choose j(r) to be nonzero only outside some sponge-free region of radius R, and when 
comparing sponge layers with other boundary treatments, we compare quantities only 
in the sponge-free region r < R. 

The sponge profile function j(r) we use is a Gaussian centred at the outer 
boundary, which we choose to place at r = 3i?, 



7(r) = 70 exp 



■3R 



(B.; 



The amplitude of the Gaussian is taken to be 70 = 1. The width a is chosen so that 
j(r) 10~ 16 (the numerical roundoff error) for r sC R, which requires a < R/3. In 
our numerical example, we take R — 41. 9M and a — 13. 3M. Hence a is considerably 
larger than the wavelength A w AM of the gravitational wave, which is required in 



order to avoid reflections from the sponge layer (cf. section 17.2.3 of [43]) . Figure B2 
shows a plot of this sponge profile. 
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